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Abstract 

We study stochastic dynamics of an ensemble of N globally coupled excitable elements. Each 
element is modeled by a FitzHugh-Nagumo oscillator and is disturbed by independent Gaussian 
noise. In simulations of the Langevin dynamics we characterize the collective behavior of the en- 
semble in terms of its mean field and show that with the increase of noise the mean field displays 
a transition from a steady equilibrium to global oscillations and then, for sufficiently large noise, 
back to another equilibrium. Diverse regimes of collective dynamics ranging from periodic sub- 
threshold oscillations to large-amplitude oscillations and chaos are observed in the course of this 
transition. In order to understand details and mechanisms of noise-induced dynamics we consider 
a thermodynamic limit ^ cxd of the ensemble, and derive the cumulant expansion describing 
temporal evolution of the mean field fluctuations. In the Gaussian approximation this allows us 
to perform the bifurcation analysis; its results are in good agreement with dynamical scenarios 
observed in the stochastic simulations of large ensembles. 

PACS numbers: 05.40-a, 05.45.Xt 
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Recent studies have shown that noise may greatly enrich dynamics of non- 
hnear systems. This contradicts an intuitive conception on a simple blurring 
effect of noise added to a dynamical system. Thus, the quantitative measures 
of noise become additional control parameters of the system, extending its pa- 
rameter space. For instance, variation of noise intensity may lead to qualitative 
transformations in the spatio-temporal dynamics. Thus noise changes the com- 
plexity of a system. As an example of such behavior we present an ensemble of 
globally coupled excitable elements. We let Gaussian white noise act additively 
and independently on each element of the ensemble. Simulations of the coupled 
Langevin equations as well as the bifurcation analysis of the equations for the 
lowest five cumulants show a rather complex sequence of qualitative changes 
when the intensity of the additive noise is varied. 



I. INTRODUCTION 



Noise-enhanced order in nonlinear dynamics far from equilibrium has become a rapidly 
growing area of modern statistical physics. Stochastic resonance directed fluxes in 
stochastic ratchets P , noise- induced phase transitions ^ are intensely discussed phenomena 
in nonlinear stochastic systems. Common to all of them is that an increase of the noise level 
results in a more ordered spatio-temporal response. This growth of order can be expressed 
and measured using various quantities. In particular, application of information theory al- 
lowed to quantify directly the noise-enhanced ordering and information transfer in bistable 
and excitable systems 14,15,0]. 

A broad class of dynamical systems with diverse applications is characterized by excitable 
dynamics. Chemical reactions js] , lasers , models of blood clotting [9] , cardiac tissues [l^ 
and nerve cells belong to its most prevalent realizations. An excitable system possesses 
a stable equilibrium state (or rest state) from which it can temporarily depart if disturbed 
by a large enough stimulus: the system responds with a spike, a large excursion in its phase 
space, returning back to the rest state through a refractory period. Fundamental phenomena 
in excitable systems, such as pulse propagation, spiral waves, spatial and temporal chaos 



and synchronization are well studied 
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Analysis of the influence of noise on excitable dynamics has recently attracted great 
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y,y,i27i as a necessary step towards realistic 



description of natural systems. The inclusion of fluctuations into models of excitabilit y is 
desirable due to their omnipresence in processes of signal transmission and detection 28] 
which constitute often an interplay of signals and noise. Hence, the strength of variability as a 



22j | . That is why 



measure of the present noise enters the excitable models as a new; parameter 
excitable dynamical systems per se are in nonequilibrium and the acting noise is unbalanced 
with respect to dissipative forces. Hence, a variation of noise often may induce qualitative 
changes in the performance and functioning of excitable systems which can be characterized 
. of ..e dynamic, of ..e. Mf.ca«on. fly 0HHQ. 

A variety of models was proposed to describe different regimes of excitable systems, 
including dynamics of excitable, spiking, and bursting neurons j35|. Here we concentrate on 
one of the most popular systems, the FitzHugh-Nagumo model which was derived from the 
Hodgkin-Huxley model of the giant squid axon [l3^. Despite its relative simplicity it displays 
a variety of biologically relevant dynamical regimes if taken as a single element or embedded 
into a network. It has been shown that the influence of noise on excitable systems such as 
the FitzHugh-Nagumo model can lead to a coherent temporal response via the phenomena 



of stochastic 1] and coherence |^ |iY|| resonance. 

In this paper we study noise- induced dynamics in a network of globally coupled FitzHugh- 
Nagumo elements, each with its own noise source. Most studies on noise-induced transitions 
co.id.ed .-called .uMpUca.ve o. pa.a.e.nc noi. o. „oi. be., colored B . In con.a., 
in our case excitable elements are driven by additive Gaussian white noise. Being statistically 
independent in different elements, noise sources possess the same intensity, T, which we use 
as a control parameter. We concentrate on the parameter region where a single uncoupled 
element possesses a stable equilibrium and is excitable. We simulate the stochastic equations 
and show that variation of T results in a rich variety of regimes of the collective response the 
ensemble, expressed in terms of the mean field (Section |H]). In order to delineate dynamical 
mechanisms behind these transitions we use cumulant approach in Gaussian approximation, 
thereby passing from a description in terms of coupled Langevin equations to a deterministic 
system of the fifth order for the cumulants of the ensemble distribution. This system is 
further analyzed by methods of the bifurcation theory fSection llII|) . In particular, separation 
of timescales allows to identify the transition from irregular minor oscillations to regime of 
intermittent spikes of large amplitude as the canard explosion of a chaotic attractor. We 
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demonstrate that the spiking regime, whether regular or irregular, is possible: (a) only in the 
restricted interval of noise intensities, and (b) for moderate values of the coupling strength: 
a too strong coupling holds the ensemble together and keeps it close to the equilibrium; a 
too strong noise maintains the uniform distribution with non-oscillating mean values. 



II. LANGEVIN APPROACH: NUMERICAL SIMULATION 

As a prototype of an ensemble of excitable units we investigate the set of noise-driven 
FitzHugh-Nagumo systems. The individual systems are coupled through the mean field and 
obey the equations: 

exi = Xi - Y -Vi + l - ^i), 



Vi = Xi + a + V2TUt), t = l,...,N. (1) 

In the biochemical context the variables Xi denote e.g. the values of the membrane po- 
tential whereas yi are responsible for the inhibitory action. Here ^i{t) are the independent 
((^j(t)^j(t)) = 6ij) white Gaussian noise sources with vanishing average which individually 
act on the recovery variables yi(t). This could be interpreted, e.g., as a random_ opening of 
ion channels 



38| which stochastically changes the conductivity of a membrane j39|. 
In the absence of noise (T = 0) the system (P) possesses a unique equilibrium state with 
Xi = —a and yi = a^/3 — a. The equilibrium is stable for > max(l, 1 — 7) and unstable 
otherwise. This follows from the characteristic equation which reads 

and means that positive coupling strength 7, whatever large it is, does not influence the 
stabihty of equihbrium. 

As instantaneous global characteristics of the set of FitzHugh-Nagumo systems we con- 
sider the first moments of the distribution: mean field components (x)(t) and {y){t). 

The evolution equations (0) were solved numerically for the ensemble of = 10^ oscil- 
lators at different values of the noise intensity T. For our computations, we take e = 0.01, 
separating thereby the characteristic timescales of the variables. Further, we choose the 
value a = 1.05 under which an individual system is excitable: if initial perturbation exceeds 
a certain threshold value, the system "fires" i.e. displays a spike of large amplitude before 



finally settling down to the steady state. Since (1 — a^)^ < 4e, the leading eigenvalues, 
which characterize the perturbation decay in the immediate vicinity of the equilibrium, are 
complex. 

We begin with the moderate value of the coupling strength 7 = 0.1 and investigate the 
temporary patterns of the mean fields for different values of the noise intensity T. 

As soon as the noisy perturbations are introduced, the individual elements, from time 
to time, are kicked across the excitation threshold and exhibit spikes. However, under 
sufficiently low values of T such irregular spikes are very seldom and do not contribute to 
the mean fields. Of course, since the noise is white, the nearly simultaneous "firing" of 
large subsets of elements should occur for arbitrarily small T, but such events appear to be 
extremely rare, and we never observed them within the integration intervals of t < 10^. 

As seen in Fig^ at low noise intensities the mean field exhibits only minor excursions 
from some average value. For (x) the time average virtually does not decline from the 
equilibrium state —a whereas the time average of (y) gets shifted. The amplitude of these 
excursions grows with the increase of T. Besides, qualitative changes occur: for very small 
T dynamics of the mean field is akin to a non-biased random walk (Fig^a)), at slightly 
larger values of T the motion gradually acquires the character of noisy rotation (Fig^b)), 
and still higher values of T yield phase portraits which are reminiscent of the spiral chaotic 
attractor (Fig^c)). In order to exclude the numerical artefacts and diminish the finite-size 
effects we performed several control runs with a much higher (A^ = 3.2 x 10^) number of 
coupled systems and obtained in the latter case the even more pronounced pattern of the 
multi-band chaotic attractor. 

Starting from T ^ 2.76 x 10~^ we observe the qualitative transition: minor ("subthresh- 
old") oscillations of the mean field are interrupted by large outbursts (spikes). This is a 
consequence of the collective phenomenon: for a large proportion of the individual units the 
firing events occur almost synchronously. Immediately beyond the threshold the spiking of 
the mean field is intermittent: the length of time intervals between the spikes is strongly 
varying (FigEK)- With the increase of T the interspike intervals become shorter and more 
homogeneous (Fig|2l3). The number of minor oscillations between two subsequent spikes 
steadily decreases (Fig|2t), proceeding from very large values (dozens in Fig|2K) to 2-3 in 
Fig|2t, and 1-2 in Fig|21i. Finally, such minor oscillations completely disappear and the mean 
field exhibits frequent non-interrupted spiking (FigI2h)- Hence variation of noise creates a 
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FIG. 1: Noise-induced onset of local oscillations around the equilibrium for a = 1.05, e = 0.01, 
7 = 0.1, N = 10^: (a) T = 10"^, (b) T = 2.4 x IQ-^, (c) T = 2.7 x 10^^. Changes in the phase 
portrait indicate a bifurcation from disordered fluctuations to subthreshold oscillations. Note the 
smallness of the oscillation amplitude. 

coherence resonance of the coupled ensemble. 

The evolution of the phase portrait for the mean field is presented in FiglHl the attractor 
which corresponds to the state of intermittent spiking consists of small and large loops 
(FigOK); the attractor of the non-interrupted spiking state is almost indistinguishable from 
the large-scale limit cycle (Figl^b)- 

Noteworthy, transition to the spiking regime is accompanied by drastic changes in the 
power spectrum of the mean field. Just below this transition the spectrum possesses a well- 
defined peak which sits upon the noisy background and corresponds to the typical frequency 
of one minor rotation around the equilibrium (Fig|3K)- Above the transition the mean field 
is characterized by the broadband spectrum (FigElD); the latter peak can still be vaguely 
identified, but there is no sharp contrast to the neighboring frequency values. A new feature 
of the power spectrum are the characteristic deep minima at the value which corresponds 
to the inverse "recovery time" as well as at its harmonics 

In fact, this minima, albeit not so pronounced, can be recognized already in Fig 13^. This 
is explained by the fact that the individual units of the ensemble exhibit intermittent spiking 
even before the transition of the mean field. Close before the threshold small clusters of such 
elements start to spike cooperatively, making thereby the inverse recovery time visible in 
the spectrum already at this stage. 

The range of noise intensities which enable the regime of spiking is relatively narrow. 
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FIG. 2: Transition from intermittent spiking of the mean field to the regular spiking pattern. The 
sequence of pictures shows an increase of coherence in the response of the ensemble with respect 
to increasing noise. The noise intensity is varied: (a) T = 2.77 x 10^^, (b) T = 2.8 x 10~^, 
(c) T = 2.9 X 10^^, (d) T = 3.0 X 10"^, (e) T = 3.1 x lO""^, other parameters like in Fig. [H 

It is convenient to characterize the oscillatory states by means of the magnitude: distance 
between the global extrema of one of the mean fields, i.e. d = max((x)) — min((x)). In 
FigEt we present the dependence d{T)\ it can be seen that after the initial abrupt growth 
the magnitude rapidly subsides. For the values of T > 0.01 no significant oscillations can 
be observed: mean fields (x) and {y) display only weak fluctuations around constant values. 
This indicates that the system has arrived at the steady distribution; however, in contrast 
to the sharp steady distribution at T = 0, in the latter case the width of the distribution 
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FIG. 3: Phase portraits of the mean field: (a) Intermittent spiking with T = 2.77 x 10 ^ , 
(b) Regular oscillations with T = 3.1 x 10"'^, other parameters like in Fig. ^ 



0) 

o 

Q. 



100 








-{ 


(a) 


1 


J' \ . 

/ v. 


0.01 









4 (0 6 



100 



1 . 



(b) 



4 (0 6 



8 



FIG. 4: Power spectra of the mean field {x){t): (a) Regime of irregular subthreshold oscillations 
at T = 2.73 X 10^^, (b) Regime of intermittent spiking at T = 2.76 x 10""^, other parameters like 
in Fig^ Minima in the spectral curve correspond to the inverse recovery time. 



is quite large. Of course, the stationary distribution does not imply stationary states for 
individual units: since the noise intensity is high, each of them remains in the regime of 
frequent spiking. However, coherence between the individual spike events is lost, and the 
ensemble average displays no time dependency. 

An increase of the coupling strength 7 produces a quantitative change in this picture 
(Fig|S|D): now the spiking state appears somewhat later and persists in a larger interval of 
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FIG. 5: Resonance-like dependence of the magnitude of oscillations on the noise intensity T for 
different coupling strength. Coupling enlarges the range of noise intensities where the resonance 
occurs. Left panel: 7=0.1; right panel: 7=1, other parameters like in Fig. ^ 

the values of T. Qualitatively we observe the similar resonance-like picture with change 
of coupling strength: there exists a certain range of noise intensities which allows for the 
spiking regime of the mean field; both the too weak noise and the too strong one suppress 
the spiking. However, the enlargement of the noise range where spiking occurs proves that 
the array enhances the collective spiking. 

The further growth of 7 leads to disappearance of the spiking regime as well as of the 
state of minor oscillations: for 7 > 2.2 we failed to observe significant temporal variations 
of the mean field at arbitrary noise intensities. Apparently, too strong coupling among the 



elements counteracts the noise: it holds individual systems together in t 
equilibrium and does not allow them to escape across the firing threshold 
1.O1 
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FIG. 6: Domains of existence of nonstationary regimes on the coupling-strength 7 vs noise intensity 
T plane: there are no spiking states outside the inner curve li and no oscillatory states at all outside 
the outer curve l2- Left panel: global view; right panel: enlarged part for small values of 7, other 
parameters like in Fig. ^ 
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Boundaries of existence of nonstationary states on the parameter plane are shown in 
Figini Here we see that not only the strong but also the weak coupling does not benefit the 
spiking state: apparently, for 7 < 0.03 the coupling is insufficient in order to synchronize 
the firing events of individual elements. It should be noted that the upper part of the curve 
li is somewhat ambiguous: here the spiking state arises from subthreshold oscillations in 
the course of decrease of T not via the abrupt transition, but rather through the continu- 
ous (although quite fast) growth of their amplitude; the plotted curve corresponds to the 
parameter values which ensure that the magnitude d equals 1. 

Summarizing the results of numerical simulations, we state that the nonstationary mean 
fields can be observed only in a restricted range of the noise intensity, provided that the 
coupling strength does not exceed a certain value. Among those nonstationary states we 
can distinguish the chaotic subthreshold oscillations, the intermittent spiking regime and 
the regular spiking pattern. 

III. GAUSSIAN APPROXIMATION: DYNAMICS OF CUMULANTS 
A. Governing equations and geometry of phase space 

The more detailed analysis of the set of coupled systems requires the deeper knowledge 
of the instantaneous distributions of Xi and In the limit of — 00 such knowledge can 
be obtained from the analysis of the corresponding Fokker-Planck equation; presence of the 
coupling term in Eq.((TI) adds to this equation the quadratic nonhnearity and destroys its 
variational character. Recently, the non- stationary Fokker-Planck equation for this problem 
was treated with the help of the expansion of the distribution density into the Hermite 
polynomials Q|; several interesting results were reported, but the global diagram of states 
is still missing. 

Description in terms of the moments of distribution turns the problem into an infinite 
set of ordinary differential equations. In order to keep the problem tractable, one should 
either truncate the set of moments or come up with a plausible closure hypothesis. We 
choose the latter way and approximate the simultaneous state of the system by the Gaussian 
distribution of Xi and yi with time-dependent parameters. Since for the Gaussian distribution 
(and only for it) all cumulants of the high order vanish, this assumption allows us to reduce 
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the infinite-dimensional problem to the low- dimensional model. 




FIG. 7: a) 3-dimensional Projection of the twodimensional slow surface in the 5-dimensional phase 
space of the cumulant equations for 7 = 0.1. Shaded domain: repehing region where inequality Q 
holds, b) Repelling region on the slow surface for 7 = 0.1, 7 = 0.5 and 7 = 1.5. 

Let us average the evolution equations over the ensemble. Conditions of vanishing of 
cumulants of the 3rd and 4th order provide expressions for the higher-order moments (x?), 
{xlui) etc. On substituting these expressions into the averaged equations, we arrive at the 
dynamical system which governs the behavior of the cumulants tjIx = (xj), rriy = (i/i), 
Dx = {{xi - ruy = {{yi - rriyY), D^y = {{xi - mx){yi - niy)) of the distribution: 
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dt 
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Dy + eDx 





For small values of e the separation of time scales in Q turns mx, Dx and Dxy into "fast" 
variables whereas my and Dy evolve on the slow time scale. Accordingly, in the 5-dimensional 
phase space there is a two-dimensional surface which corresponds to slow motions. 

Location of this surface in the limit of vanishing e is obtained by setting e = in Eq.Q. 
It is convenient to parameterize the surface by coordinates mx-, Dx'- given arbitrary values 
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of these two variables (of course, only non- negative values of are physically meaningful) 
the three remaining coordinates are determined as 

my = rrin, - ml/?> - mn^Dx 

D,y = Dx{l - D., - ml - ^) (3) 

2 



Dy = D,{l-D.,-mi--f) 



Solving the linear system 



[1 — Dx — ml) rfix —m^Dx = m^ + a 

T 

-(4m^) Dx rhx +(1 - -ml- -f)Dx = Dx + (4) 

1 - D^. - - 7 

for variables nix and Dx, yields evolution equations for dynamics upon the slow surface. 
[Explicit expressions are straightforward, but too long to be quoted here]. Stability of this 
surface with respect to transversal perturbations is governed by the sign of the determinant 
of the system (@)): slow surface is repelling for 

{ml - ly + DxiSDx - 4) < 7(1 -ml- Dx) (5) 

and attracting otherwise; notably, position and shape of the repelling region depend only 
on one of the parameters: the coupling strength 7. 

One of the projections of the slow surface in the phase space is shown in FigjT^; the 
butterfly-shaped repelling region is shaded. The geometry of the phase space has impli- 
cations for the dynamics in general. We can expect that phase trajectories approach the 
attracting part of the slow surface and move along it; in case of penetrating into the region 
dni), they move for a certain time along the repelling surface (effect known as the "canard 
phenomenon"), depart from the slow surface and eventually approach the attracting part 
again. 

Indeed, the bifurcation analysis of Eq.(j21) discloses this picture. 
B. From stationary distribution to small-scale chaos 

We consider the equations (j21) in the parameter range a > 1 where the solitary FitzHugh- 
Nagumo oscillator has a stable steady state. For T > the equations (j2)) also possess the 
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unique steady solution in the domain Dx,Dy > 0: 




my = — -a + aD^, Dy = eD^ + T{a^ + + - 1), (6) 

O 

D,y = -T 

(there exists also the physically meaningless unstable state with negative and Dy). This 
equilibrium describes the stationary distribution of oscillators; in the phase space the cor- 
responding fixed point is located on the slow surface. 

Destabilization of the equilibrium © under the increase of the noise intensity has a 
convenient graphical interpretation. Recall that shape and position of the repelling region 
on the slow surface depend only on 7. FigI7|D shows the outline of this region for several 
values of 7. As seen from Eq.®, the equilibrium value of rrix is T- independent, whereas 
the value of is the monotonically growing function of T. At T=0 the fixed point lies on 
the abscissa of FigEb- For all values of 7, the border of the repelling region passes through 
the point m^.=-l, Dy=0] part of the abscissa to the left from this point always lies in the 
attracting region. Accordingly, for a > 1 and arbitrary 7 the stationary solution is stable in 
the absence of noise and (by continuity) at very small values of T. 



stable 
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FIG. 8: (a) State diagram of the equilibrium at a = 1.05, e = 0.01. Region of stability is bounded 
by the curve of the Andronov-Hopf bifurcation, (b) Bifurcation surface of the Andronov-Hopf bi- 
furcation in the parameter space of Eqs© with e = 0.01. The shaded region indicates a subcritical 
bifurcation. 

As a consequence of the increase of T, the fixed point moves upwards along the vertical 
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line TJix = —a and can cross the repelling region. Destabilization occurs when the fixed point 
hits the lower border of this region; in terms of Eqs Q this event is the Andronov-Hopf 
bifurcation. In the limit e — the critical value Th is given by the smaller root of the 
quadratic equation 



where h = a? — 1. 

To enable this instability, the vertical line m^. = —a should cross the repelling region in 
FigEb. In terms of the parameter values this condition reads as 



For 7 > 7o the fixed point on the slow surface is always located to the left of the repelling 
region, and the equilibrium is stable for all values of T. This can be interpreted as stabi- 
lization of the equilibrium distribution by strong coupling. With the increase of a the value 



of 7o becomes smaller; it vanishes at a = = y 1 + y4/3 ~ 1.467. For a > the steady 
state is stable indendently of the values of 7 and T. At small nonzero e the qualitative and, 
largely, the quantitative picture persists; merely the values of 70 and slightly change. 

Under sufficiently high values of T the fixed point moves over the upper border of the 
repelling region on the slow surface: the equilibrium regains stability. Again, this is the 
Andronov-Hopf bifurcation; the asymptotics at e ^ for the corresponding value of T is 
given by the larger root of Eq.(|7j). Shape of the bifurcation curve on the plane of parameters 
7 and T is shown in FiglHK- 

Thus we observe that in the domain of low and moderate values of 7 the equilibrium 
exhibits two Andronov-Hopf bifurcations: the destabilizing and the stabilizing one. Compu- 
tation of the cubic term in the amplitude expansion discloses a relatively small area on the 
parameter plane of a and 7 where the first bifurcation is subcritical (Fig|H|D): the unstable 
periodic solution branches off in the direction of small T. Here, a hysteresis occurs: the 
time-independent regime coexists with the stable oscillatory state; the latter is born from 
the tangent bifurcation shortly before the destabilization of the equilibrium. For most of 
the values of a and 7, however, the first Andronov-Hopf bifurcation is supercritical: the 
stable limit cycle with small amplitude ~ \/T — Th and frequency ~ e^^^^ is born. The 
second bifurcation is always supercritical: the stable limit cycle is born into the domain of 



9T^ + T{16b^ - 16 - 126 - 47 + 967 + 27^) + 2b{b + 7)^(2 + 26 + 7) = 



(7) 



7 < 70 = 2(3a2 - 1 - 2aV3a2-3). 
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FIG. 9: Projections of phase portraits for EqsQ (cf. Fig. El for different regimes of local sub- 
threshold oscillations. The dotted grid SI shows the slow surface in the phase space. The surface 
SI is stable to the left of the dashed line B. (a) T = 0.00157, (b) T = 0.00158, (c) T = 0.0015826, 
(d) T = 0.001585. Other parameters: a = 1.05, 7 = 0.1, e = 0.01 

lower values of T. With regards to finite-amplitude states, there seem to be no non-decaying 
solutions in the part of the parameter plane above the bifurcation curve in Fig|H^. 

Near the lower branch of the bifurcation curve, the newborn limit cycle includes the 
segment (denoted by a in FigEti') which lies close to the attracting region of the slow surface, 
and the segment r which leads along the repelling region; from there, the orbit leaps back 
to the attracting region (segment / in FigEt-)- 

Since the slow surface is two-dimensional, one can expect new qualitative effects in com- 
parison to the single FitzHugh-Nagumo system with its S-shaped one-dimensional slow 
manifold. Indeed, when T is further increased, the periodic orbit loses stability by means of 
the period-doubling bifurcation. This event is followed by the sequence of further period- 
doublings which culminates in the onset of the chaotic state. In terms of the ensemble, 
the latter describes minor (localized) chaotic oscillations around the stationary distribution. 
Transformation of the attracting orbit and its position near the slow surface is depicted in 
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Figia 




FIG. 10: Transition in the power spectrum of the cumulant equations, a = 1.05, 7 = 0.1, e = 0.01. 
(a) r=0.001585, (b) T=0.001586. 

The power spectrum of a trajectory on the chaotic attractor consists of the peak at 
the characteristic frequency of the minor oscillations, its harmonics and the broad noisy 
background fFiglTUk). 



C. Prom chaotic to regular spiking 

Further growth of T brings about dramatic changes in the shape and size of the attractor. 
As mentioned above, attracting orbits involve segments which lie close to the repelling part 
of the slow surface. In such situation a minute variation of the parameter may be a reason for 
the abrupt (albeit continuous) increase of the attractor size by several orders of magnitude: 
the so-called "canard explosion" . This phenomenon has been well documented for the case 
when the attractor is a limit cycle (e.g. in the forced Van der Pol equation 42j, and in the 
context of a single FitzHugh system disturbed by noise The peculiarity of our 



situation is that the "canard explosion" happens not to an individual periodic orbit but to 
the chaotic attractor as a whole (of course, individual unstable periodic orbits embedded 
into the attractor also experience the explosion, but this remains unnoticed for an observer 
who watches only the stable sets). 
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FIG. 11: Regime of irregular spiking in Eqs ((21); a=1.05, 7=0.1, e=0.01, T=0. 001586. (a) projection 
of phase portrait; (b) enlarged part of phase portrait; (c) plot of mx{t); (d) plot of Dx{t). 

Projections of phase portraits and temporal evolution of the mean fields rux and rriy in 
this state are shown in FiglTTl Apparently, this regime corresponds to intermittent chaotic 
spiking: after many minor oscillations in the vicinity of the unstable steady equilibrium, the 
system exhibits a spike of large amplitude, followed by the next epoch of localized chaotic 
oscillations. 

A joint graphical presentation of the attractor and the slow surface is helpful for under- 
standing the nature of spiking. The enlargement of the domain around the equilibrium in 
FiglT^ discloses that each arrival of a trajectory to this domain occurs via a steep descent 
into the "valley" along the slow surface. This return is followed by several minor revolutions, 
which share the familiar pattern: two segments along the, respectively, attracting and the 
repelling regions of the slow surface, and a return back to the attracting region. Finally, the 
trajectory "gets through" the slow surface 0] and is repelled far away from the equilibrium. 
A global picture fFigJ12b) indicates that this flight ends in the domain of positive m^, where 
the trajectory rapidly approaches the slow surface and moves along it for a certain time; the 
spike is accomplished by another fast flight back to the initial domain. 

Accordingly, each single spike of the mean field is accompanied by two spikes of the 
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FIG. 12: Projection of the attractor and the slow surface in the state of irregular spiking. Left: 
global view of the attractor; right: enlargement of the domain near the unstable equilibrium. 

variance (Fig. IT^ : the latter remains small on each of two slow segments but rapidly peaks 
during each of two fast flights. 

Transition to the spiking regime is accompanied by sharp changes in the power spectrum 
of the process; this is shown in FigJlUb. 

States with chaotic spiking occupy only a narrow strip in the parameter space. Further 
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FIG. 13: Time evolution of the mean field and its variance within a single spike, (a) variable m^(t); 
(b) variable Dx{t). Parameter values like in Fig llll Rapid changes oi nix in either direction induce 
peaks of Dx- 




0.75 
X 0.5 

Q 

0.25 



18 




FIG. 14: Transition from intermittent to regular spiking in the cumulant equations with e = 0.01, 
a = 1.05, 7 = 0.1: (a) T = 0.00168, (b) T = 0.00172, (c) T = 0.00185, (d) T = 0.0018569, (e) 
T = 0.00220, (f) T = 0.0023105, (g) T = 0.0024. 

increase of T regularizes the spiking pattern: individual spikes are now separated by equal 
amounts of minor oscillations. With growth of T this amount gets smaller: it is a kind of 
the inverse "period-adding" sequence in which localized oscillations between two subse- 
quent spikes are replaced by A^-1 ones (there is also a narrow transition region with more 
complicated states), then A^-1 yields to etc Finally the minor oscillations com- 
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pletely disappear and the system proceeds into the regime of uninterrupted periodic spiking. 
Different stages of this process are shown in FiglT^ whereas the basic transition curves on 
the bifurcation diagram are plotted in FiglT^. 
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FIG. 15: a) Part of the bifurcation diagram for cumulant equations at a = 1.05, e = 0.01 for small 
values of the coupling 7. AH: Andronov-Hopf bifurcation; Canard explosion: onset of intermittent 
spiking oscillations, I3: onset of regime with 2 minor oscillations between spikes, l2'- onset of regime 
with 1 minor oscillation between spikes, h: onset of regime without minor oscillation between spikes 
b) Important transitions on the parameter plane: AH: Andronov-Hopf bifurcation; pdi 2: period- 
doubling bifurcations. 




D. Strong noise: back from regular spiking to stationary state 

As discussed above, under high values of T the only attractor of the system is the sta- 
tionary state. Formally, it is the same solution (jHl) as in the case of low T; however, now 
this state describes the broad stationary distribution with high values of and Dy. Here 
we briefly outline the way to this stationary state from the regular spiking pattern. 

Quantitatively, we characterize the oscillations in terms of the same "magnitude" as in 
Sectjn] now this "attractor diameter" is defined as d = max(m^) — min{mx). Dependence 
of d on T for several fixed values of 7 is shown in FiglTBl 

We observe that the weak decay of the amplitude at low values of T is interrupted by the 
sharp decrease by a factor of ~ 2 which occurs within the narrow interval of T. This decrease 
is followed by the final interval of continuous decay until the periodic state disappears in 
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FIG. 16: Resonance like dependence of the magnitude of oscillations as function of the noise 
intensity T. Results of integration of the cumulant equations. 

the Andronov-Hopf bifurcation; here, the segments of the curves obey the square-root law. 
The sharp decrease reminds of the canard explosion at low values of T, although it is less 
dramatic. At a closer look we see that the decrease is always preceded by a short interval 
in which d{T) is non-monotonic. Checking the corresponding parameter values, we disclose 
here additional bifurcations. The local minimum of d{T) corresponds to the period-doubling 
which occurs on the upper branch of the curve pd2 on the parameter plane in FigllSb. 

This event is followed by the whole period-doubling scenario and the narrow range of T 
in which the spiking pattern turns chaotic again. In contrast to the previously described 
intermittent chaotic spiking, there are no small oscillations between the subsequent spikes; 
only the amplitudes of the spikes differ. The magnitude d weakly grows again, but in the 
course of the further minor increase of T the loops of the orbit "slide off" the broad part of 
the repelling region of the slow surface, and d rapidly decreases. Meanwhile, the chaotic state 
turns back into the regular one, and finally the reverse period-doubling bifurcation restores 
the simple periodic state. On the parameter plane the curve of the latter bifurcation turns 
out to be the upper branch of the already familiar period-doubling curve described in the 
previous subsection (line pdi in Fig. [T3b). Fast "halving" of the attractor diameter takes 
place just below the curve pdi\ in the broad parameter range above this curve and below 
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the upper branch of the curve AH of the Andronov-Hopf bifurcation the amphtude of the 
hmit cycle slowly decays. 

For the values of 7 to the left of the curve pd2 (1.7 < 7 < 2.2) the ordering of states is 
slightly different. Here the reverse period-adding sequence which simplifies the intermittent 
spiking pattern at low values of T, ends up in a limit cycle with several (4-6) spikes which 
exists over the broad range of T. Slightly below the upper branch of the curve pdi the tran- 
sition to chaos occurs; it is followed by the reverse sequence of period-doubling bifurcations 
and very strong decrease of the diameter d. 

Summarizing the results obtained from the analysis of the equations which govern the 
dynamics of the cumulants, we observe their remarkable qualitative correspondence with the 
basic knowledge obtained from the integration of the Langcvin equations. Both the very 
weak and the very strong noise correspond to stationary probability distributions. Descrip- 
tion in terms of the cumulants also confirms the stabilization of the stationary distribution 
by the strong coupling. The attr actors of the cumulant equations practically do not differ 
from the phase portraits of the mean fields computed from simulations of the Langevin 
equations; the patterns of power spectra, especially for the state of intermittent spiking, 
also look strikingly similar. Besides, the analysis of the cumulant equations discloses the 
detailed picture of transitions between different non-stationary states: birth of the intermit- 
tent chaotic spiking from chaotic subthreshold oscillations and the subsequent regularization 
of the spiking pattern under the action of noise. Quantitatively, in the cumulant equations 
the critical values of the noise intensity which correspond to transitions between different 
states, turn out to be noticeably higher than those in the numerical simulations of large 
ensembles of FitzHugh-Nagumo systems. 

IV. DISCUSSION AND CONCLUSIONS 

We have investigated an ensemble of globally coupled FitzHugh-Nagumo elements under 
the action of additive Gaussian white noise. With the help of numerical simulations of 
stochastic equations we have verified that noise intensity becomes an additional control 

parameter of the system. By changing noise intensity we provoke qualitative transformations 
of the global response of the system, characterized by the mean field. 

In order to gain better understanding of dynamical mechanisms of these transitions we 
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employed an approximate description in terms of cumulants of the Gaussian distribution. On 
the whole, modeling a large ensemble of globally coupled excitable units by the Gaussian 
distribution with time-dependent characteristics proved to be a fruitful approach: it has 
reduced the problem to the low- dimensional deterministic dynamical system which, in the 
large part of the parameter space, reproduced basic properties of the underlying stochastic 
ensemble. This allows us to conjecture that the fine details of the transitions between the 
states, accessible only in the framework of this model, keep relevance also for the whole 
ensemble. 

In our example only the slow variables have been subjected to the additive noise. Of 
course, in the natural environment and in the laboratory experiments one cannot reduce the 
influence of fluctuations to only the slow variables: the fast ones can be affected as well. This 
can be accounted for, if the additive sources of Gaussian noise are included into the first of 
Eq.(^. We performed a set of numerical experiments with such equations, and observed, in 
general, the states similar to those described in the previous sections: steady distributions 
at high or very low noise intensities, subthreshold oscillations and spiking states in the 
intermediate range. The analysis of the respective cumulant equations can be pursued along 
the same lines as above. Since the noise acts on the fast variables, the shape and position 
of the repelling part of the slow surface depend now on the noise intensity; qualitatively, 
however, the only new effect in comparison with those described in Sect JIIII seems to be 
the direct transition from periodic subthreshold oscillations to regular spiking without the 
mediation of chaotic states; this is, of course, the usual canard explosion. 

Of course, the Gaussian distribution is merely an approximation, and its applicability has 
certain restrictions. Thus, we failed to obtain the quantitative correspondence for small val- 
ues of the coupling strength 7. Stochastic simulations indicate the absence of nonstationary 
regimes in this parameter range (cf. FigEjD) whereas the analysis of the cumulant equations 
predicts the whole bifurcation scenario for arbitrarily low values of 7 and even at 7 = (cf. 
FiglTHk). The higher the value of 7, the better the correspondence between the behavior of 
the stochastic equation and the deterministic model However, since the considered 
coupling is linear, an increase of its strength reduces the complexity of dynamics towards 
a linear system. Hence, under the strong coupling the Gaussian approximation becomes a 
more adequate model of the process, but the process itself gets simpler: there is little place 
for nonlinear events and, hence, variation of T fails to invoke transitions in the behavior of 
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the mean field. 
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